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Abstract 

We  propose  the  use  of  proper  orthogonal  decomposition  (POD)  techniques  as  a  reduced 
basis  method  for  computation  of  feedback  controls  and  compensators  in  a  high  pressure 
chemical  vapor  deposition  (HPCVD)  reactor.  In  this  paper,  we  present  a  proof-of-concept 
computational  implementation  of  this  method  with  a  simplified  growth  example  for  III-V 
layers  in  which  we  implement  Dirichlet  boundary  control  of  a  dilute  Group  III  reactant 
transported  by  convection  and  diffusion  to  an  absorbing  substrate  with  no  reactions.  We 
implement  the  model-based  feedback  control  using  a  reduced  order  state  estimator  based 
on  observations  of  the  flux  of  reactant  at  the  substrate  center.  This  is  precisely  the  type 
of  measurements  available  with  current  sensing  technology.  We  demonstrate  that  the 
reduced  order  state  estimator  or  compensator  system  is  capable  of  substantial  control 
authority  when  applied  to  the  full  system.  In  principle,  these  ideas  can  be  extended  to 
more  general  HPCVD  control  situations  by  including  multiple  species  with  gas  phase 
reactions  and  surface  reactions. 
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1  Introduction 


Production  of  advanced  optoelectronic  integrated  circuits  requires  stringent  control  of 
layer  thickness  and  composition.  These  requirements  can  be  addressed,  in  part,  through 
open- loop  optimization  [6,  21,  34],  However,  because  of  process  variability  and  the  in¬ 
creasing  demands  put  on  control  of  layer  thickness  and  composition  by  state-of-the  art 
devices,  real-time  control  of  film  growth  is  desirable  [17,  37,  35,  7]. 

Many  materials  utilized  in  today’s  electronics  industry  are  manufactured  using  chemical 
vapor  deposition  (CVD)  processes  operating  at  low  pressure.  However,  there  are  also 
materials  of  potential  industrial  use  that  can  not  be  adequately  produced  at  desirable 
process  temperatures  under  low  pressure  conditions,  e.g.,  InN  films  that  exhibit  relatively 
high  decomposition  pressure  as  compared  to  other  III-V  compounds.  In  such  cases,  it 
may  be  desirable  to  extend  the  CVD  processing  to  higher  pressures.  We  are  presently 
involved  in  a  collaboration  with  material  scientists  at  N.C.  State  University  to  design 
and  build  such  a  HPCVD  reactor  with  real-time  sensing  and  control  as  an  innovative 
feature  of  this  prototype  reactor. 

There  are  many  technical  issues  associated  with  real-time  control  of  a  HPCVD  reactor, 
such  as  coupling  of  species  transport,  gas  phase  reactions,  and  surface  reactions  to  film 
growth  and  development  of  suitable  sensing  technology.  We  propose  using  accurate 
model  simulations  of  species  transport  and  reactions  to  provide  data  about  the  state 
of  the  system  and  to  design  state  feedback  controllers.  In  general,  a  full  mathematical 
model  describing  transport  processes  in  CVD  systems  is  given  by  a  system  of  nonlinear 
partial  differential  equations  representing  the  continuity,  momentum,  energy,  and  species 
equations  of  state.  Therefore,  numerical  simulations  and  control  designs  of  such  systems 
using  finite  element,  finite  difference,  or  spectral  methods  will  lead  to  a  very  large  system 
of  ordinary  differential  equations  rendering  real-time  full  model-based  feedback  control 
design  infeasible.  To  reduce  the  dimensionality  of  the  system,  we  propose  using  the 
method  of  proper  orthogonal  decomposition  (POD)  to  more  efficiently  represent  the 
system  data  [27].  Recognizing  that  observation  of  the  full  state  of  the  system  can  not 
be  experimentally  realized  with  current  sensing  technology,  we  propose  using  a  state 
estimator  or  compensator  based  on  P-polarized  reflectance  spectroscopy  (PRS)  sensing 
technology,  which  provides  the  capability  for  real-time  observation  of  the  growth  rate  [5]. 

In  this  paper,  we  present  a  proof-of-concept  implementation  of  this  method  with  a  sim¬ 
plified  growth  example  for  III-V  layers,  which  can,  in  general,  be  extended  to  more 
experimentally  relevant  cases.  Normally,  during  the  growth  of  III-V  films  the  Group 
III  reactants  determine  the  growth  rate  and  composition,  since  the  Group  V  reactants 
are  supplied  in  overabundance,  while  the  Group  III  reactants,  possessing  high  sticking 
coefficients,  are  supplied  in  dilute  amounts  relative  to  the  carrier  gas.  We  present  im- 
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plementation  of  Dirichlet  boundary  control  of  a  dilute  Group  III  reactant  transported 
by  convection  and  diffusion  to  an  absorbing  substrate  with  no  reactions.  Computational 
fluid  dynamics  (CFD)  simulations  provide  data  on  the  full  system  and  are  used  to  con¬ 
struct  a  POD  reduced  order  model. 

The  control/compensator  problem  (or  LQG  tracking  problem)  is  formulated  as  a  linear 
quadratic  regulator  (LQR)  tracking  problem,  where  the  state  of  the  system  is  estimated 
via  a  compensator  gain  from  observation  of  the  flux  of  the  reactant  to  the  substrate 
(assumed  proportional  to  the  growth  rate),  with  the  observed  flux  attempting  to  track  a 
desired  flux  value.  Dirichlet  boundary  control  at  the  inlet  determines  the  mass  fraction  of 
the  incoming  Group  III  reactant.  We  demonstrate  that  the  reduced  order  state  estimator 
or  compensator  system  is  capable  of  substantial  control  authority  when  applied  to  the 
full  system.  In  principle,  these  ideas  can  be  extended  to  more  general  HPCVD  control 
situations  by  including  multiple  species  with  gas  phase  reactions  and  surface  reactions. 

The  utility  of  POD  (also  known  under  the  names  of  principal  component  analysis  [19] 
and  Karhunen-Loeve  expansion  [20,  24])  as  a  method  of  feature  extraction  is  well  known 
in  statistical  and  pattern  recognition  fields  [16]  and  has  been  applied  in  such  diverse  areas 
as  materials  processing  [27]  and  characterization  of  human  faces  [30].  The  POD  method 
is  a  linear  transformation  of  a  multivariate  data  set  into  an  optimal  set  of  uncorrelated 
variables  (POD  modes).  The  original  multivariate  data  can  be  written  as  linear  combi¬ 
nations  of  the  POD  modes.  In  many  cases  the  POD  modes  more  efficiently  describe  the 
variability  of  the  original  data  and  some  dimensional  reduction  is  possible  by  retaining 
only  the  most  important  modes. 

POD  based  approximation  methods  have  been  widely  discussed  in  the  context  of  turbu¬ 
lent  coherent  flows  [3,  12,  13,  15,  18,  25,  31]  -  see  also  the  surveys  [11,  26].  More  recently, 
the  possibility  of  POD  based  control  design  has  been  proposed.  In  particular,  the  first 
uses  of  POD  approximations  in  optimization  or  open  loop  control  were  demonstrated  in 
[27,  34]  for  model-based  methods  and  subsequently  in  [28]  for  model-free  methods.  The 
first  use  of  POD  approximations  in  feedback  control  design  was  reported  in  [8,  9]  for 
control  of  periodically  disturbed  structures.  POD  methods  for  state  feedback  were  also 
used  at  approximately  the  same  time  in  [22]  in  an  example  for  Burgers’  equation  and 
were  later  used  in  [1]  in  control  design  for  the  heat  equation.  To  our  knowledge,  the 
results  reported  below  offer  the  first  computational  evidence  of  the  successful  use  of  POD 
based  methods  for  feedback  control  and  compensator  design  in  tracking  problems  for  a 
flow  related  system  where  only  partial  state  observations  are  available:  in  this  case  in 
the  control  of  pulsed  metallic  vapors  in  a  CVD  reactor. 

There  are  a  number  of  nontrivial  issues  related  to  the  use  of  reduced  basis  methods  in 
general  and  POD  methods  in  particular  as  a  foundation  for  approximation  methods  in 
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infinite  dimensional  systems  such  as  those  modeled  by  distributed  parameter  or  partial 
differential  equation  systems.  The  foremost  revolves  around  whether  the  infinite  dimen¬ 
sional  system  itself  can  be  approximated  well  by  a  finite  span  of  basis  elements.  That 
is,  are  the  features  of  importance  in  a  given  investigation  of  the  system  essentially  finite 
dimensional  in  nature?  There  is  growing  evidence  that  the  answer  to  this  question  is 
positive  for  many  structural  systems  and  for  a  substantial  number  of  fluid  and  electro¬ 
magnetic  applications.  Even  when  this  question  can  be  answered  in  the  affirmative,  it 
is  not  at  all  clear  that  one  can  use  such  an  approach  for  control  design  or  in  inverse 
problems  related  to  damage  detection  (another  area  we  are  currently  investigating  with 
POD  based  methods). 

In  particular,  if  one  uses  ‘snapshots’  of  the  uncontrolled  system  (as  was  done  in  [8,  9]) 
to  construct  the  POD  basis  elements,  there  is  little  reason  to  expect  that  control  design 
based  on  this  finite  dimensional  approximation  will  be  effective  when  applied  to  the 
original  system.  The  controlled  original  system  itself  may  require  a  different  set  of  finite 
dimensional  elements  for  efficient  approximation  (or,  worse  yet  may  not  be  amenable  to 
a  low  order  basis  approximation).  The  results  of  [8,  9]  suggest  that,  at  least  in  some 
situations,  happily  this  is  not  the  case.  The  approach  we  offer  in  this  paper  (which  was 
also  successfully  used  in  open  loop  control  problems  in  [27])  illustrates  that  another  tactic 
can  be  effective.  That  is,  one  may  wish  to  ‘snapshot’  on  the  system  under  several  levels 
of  control  input  (not,  in  general,  derived  from  any  optimal  or  suboptimal  design)  in  order 
to  take  snapshots  on  the  system  under  nontrivial  control  inputs. 

In  Section  2  below  we  describe  the  particular  system  under  investigation  here  and  describe 
the  simulations  needed  to  develop  POD  based  control  design.  A  sample  of  some  of  our 
computational  results  based  on  this  approach  is  given  in  Section  3.  Brief  conclusions  are 
then  summarized. 


2  Methods 

2.1  Governing  Equations 

The  governing  equations  for  the  transport  dynamics  in  CVD  processes  are  described  by 
the  following  system  of  partial  differential  equations: 

(mass) 

|Tv.(H  =  0.  (1) 
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(momentum) 

P  +^- Vnj  = -VP  +  V  •  f  -  (2) 

where  the  viscous  stress  tensor  is  of  the  form 

r  =  -^//(V  •  v)I  +  ju(W  +  WT),  (3) 

o 

(energy) 

f>cp(^+?'Vr)  =V.(tVT),  (4) 

and  (species) 

dY  1 

—  +  n  •  VY  =  -V  •  (pDVY )  ,  (5) 

where  in  (5)  we  assume  that  there  are  no  reactions  in  the  gas  phase.  In  addition,  g  is 
the  gravitational  acceleration,  v.  T,  and  P  are  the  velocity,  temperature,  and  pressure, 
/i,  cp,  and  k  are  the  viscosity,  specific  heat,  and  conductivity  of  the  carrier  gas,  D  is 
the  diffusivity  of  the  reactant,  and  Y  is  the  mass  fraction  of  the  reactant.  The  density 
variations  are  modeled  as  p  =  po[l  —  P(T  —  To)],  where  To  is  a  reference  temperature,  po 
is  a  reference  density  calculated  from  the  ideal  gas  law  at  the  reference  temperature  and 
reactor  pressure,  and  f3  is  the  volume  coefficient  of  expansion  (/3  =  1/T).  The  boundary 
conditions  for  the  above  system  of  equations  (l)-(5)  will  be  given  in  subsequent  sections. 

We  summarize  our  use  below  of  the  system  (l)-(5)  in  numerically  investigating  reduced 
order  control  and  compensator  design.  Since  only  a  trace  amount  of  Group  III  reactant 
is  mixed  with  the  carrier  gas,  we  can,  first,  solve  (l)-(4)  for  steady  state  solutions  using 
a  commercially  available  CFD  finite  element  software  package.  These  values  are  then 
employed  in  a  transient  simulation  of  (5)  to  solve  (again  using  standard  finite  element 
procedures)  for  values  to  be  used  as  ‘simulated  experimental  data’  or  ‘time  snapshots’ 
for  (5).  We  then  carry  out  the  POD  procedure  described  below  to  produce  POD  basis 
elements  {H h}iLi  using  POD  nodal  values  {z^}^=  1.  The  control  problem  for  (5)  is  then 
approximated  in  penalty  form  (10),  which  is  then  written  in  weak  or  variational  form 
(11).  A  standard  Galerkin  (quadratic)  finite  element  method  is  used  to  approximate 
solutions  of  (11);  this  yields  the  full  dimensional  approximation  of  (11)  that  plays  the 
role  of  a  reactor  simulator.  A  reduced  order  (M  <<  K)  POD  basis  is  used  to  approximate 
solutions  of  (11)  via  a  Galerkin  procedure  which  yields  the  finite  dimensional  (low  order) 
approximation  of  (11)  to  be  used  in  LQG  control/compensator  gain  design. 
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2.2  CFD  simulations 


Specifically,  we  use  a  2D  rectangular  geometry  (Fig.  1)  representing  the  longitudinal 
cross  section  through  the  center  of  an  HPCVD  reactor  with  dimensions  (height=0.011 
m,  length=0.156  m,  substrate  length=0.048  m)  that  are  similar  to  the  present  N.C.  State 
reactor  design. 


n 


n 


r6 


SUBSTRATE  15 


r4 


r3 


X 


Figure  1:  HPCVD  reactor  geometry. 

The  governing  equations  are  discretized  using  the  Galerkin  finite  element  method  with 
weighted  residuals  for  the  degrees  of  freedom  (v.  P,  and  T ).  A  mixed  formulation  with 
132  quadrilateral  elements  (corresponding  to  453  nodes)  is  used  with  piecewise  linear  dis¬ 
continuous  elements  for  pressure,  and  quadratic  (8-noded)  elements  for  the  other  degrees 
of  freedom.  Solutions  for  the  variables  v,  P,  T,  and  Y  are  obtained  using  simulations 
with  commercially  available  code  (FIDAP,  Fluid  Dynamics  International,  Evanston,  IL) 
on  a  Silicon  Graphics  Power  Indigo  2. 

Steady-state  simulation.  Under  the  dilute  approximation,  steady-state  solutions  for 
the  velocity,  temperature,  and  pressure  are  independent  of  the  reactant  concentration 
and  are,  of  course,  time- independent.  Therefore,  the  spatially-dependent  solutions  for 
v  and  T,  appearing  in  (5),  are  obtained  from  a  steady-state  simulation  of  (l)-(4)  at 
atmospheric  pressure  with  hydrogen  carrier  gas.  Temperature  dependent  values  for  /i,  k , 
and  cp  necessary  for  solution  of  (l)-(4)  are  linearly  interpolated  from  measurements  taken 
from  the  available  literature  [23,  32,  36].  A  parabolic  velocity  flow  profile  is  specified  at 
the  inlet  (n),  with  an  average  inlet  velocity  of  0.1147  m/s.  Room  temperature  boundary 
conditions  are  imposed  at  the  inlet  and  along  the  upper  wall  (T2).  Along  the  bottom 
wall,  the  substrate  (r5)  temperature  is  fixed  at  800°K  ,  with  a  non-linear  temperature 
decrease  from  the  substrate  edge  to  the  inlet  (T6)  and,  similarly,  from  the  substrate  edge 
to  the  outlet  (r4). 

Reactant  transport  simulation.  Trimethylindium  is  an  organometallic  precursor  com¬ 
monly  used  as  a  source  material  for  CVD  growth  of  III-V  films  where  indium  is  one  of 
the  constituents  ( e.g .,  InN).  A  transient  simulation  of  trimethylindium  (TMI)  convection 
and  diffusion  (5)  is  used  to  obtain  data  for  construction  of  the  POD  modes.  This  system, 
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including  a  nonzero  boundary  input,  is  given  by 


dY 

~dt 


+  v-VY 
Y (0,  x) 
Y  (t,  x) 
Y (t,  x) 
dY(t,  x) 
~&n 


iv  .  ( pDVY ) 

0 

1  on  n 

0  on  T5 

0  on  T2  U  T3  U  T4  U  T6, 


(6) 


pj 

where  Y  is  the  mass  fraction  of  TMI  and  -Y-  denotes  the  outward  normal  derivative. 

on 

Initially,  the  TMI  mass  fraction  is  assumed  zero  everywhere  in  the  reactor.  A  normalized 
TMI  mass  fraction  of  one  (Y  =  1)  is  maintained  as  a  nontrivial  input  ‘control’  at  the  inlet 
(ri).  The  reactor  walls  (T2,  F4,  and  T6)  are  non-absorbing,  and  the  substrate  (T5)  is 
assumed  to  be  perfectly  absorbing  (concentration  of  zero).  Values  for  v  and  p  appearing 
in  (5)  are  provided  by  the  steady-state  solutions.  Temperature  dependent  values  for  the 
diffusivity  D  of  TMI  in  hydrogen  are  linearly  interpolated  from  values  taken  from  the 
available  literature  [33].  Transport  of  the  TMI  is  simulated  using  standard  finite  elements 
over  a  2-s  time  period,  sufficiently  long  for  achieving  a  steady-state  concentration  in  the 
reactor.  Time  integration  is  implemented  using  a  backward  Euler  method  with  fixed  time 
steps  (0.02  s).  Intermediate  solution  values  are  stored  at  each  time  step  to  be  used  later 
in  construction  of  the  POD  basis  elements. 


2.3  Construction  of  POD  Modes 

The  reactant  transport  simulation  described  above  provides  a  multivariate  data  set  con¬ 
sisting  of  K  vectors  X  =  {x^,  x£, . . . ,  x ^},  each  vector  representing  N  nodal  values  of  the 
species  concentration  in  the  reactor  at  different  times  during  the  period  that  the  reactant 
was  entering  the  reactor  and  reaching  its  equilibrium  distribution.  In  this  calculation, 
100  time  snapshots  are  available  but  only  iv  =  60  are  used  in  the  POD  construction  (see 
Figure  3  and  the  associated  discussion  below).  This  original  data  set  X  is  transformed 
to  a  new  set  of  uncorrelated  variables  (POD  modes) 

Z  =  {z?,z?1...,z$}  =  X*,  (7) 

where  the  columns  of  $  =  {4>f,  4>2  ,  ■  ■  ■ ,  4>k}  are  the  eigenvectors  of  the  product  matrix 
(X'X =  A i<j)f ,  ranked,  in  descending  order,  with  respect  to  the  associated  eigen¬ 
value.  The  prime  superscript  denotes  the  transpose  of  the  matrix.  The  POD  modes  Z 
are  orthogonal  zf  ■  z ^  =  XiSij,  and  the  transformation  of  variables  preserves  the  data 
variability 

K  K  K 

£(*'*)«  =  J2(z'z)kk  =  E A*-  (8) 

k  k  k 
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Expansion  of  the  original  data  X  in  terms  of  the  most  significant  POD  modes  minimizes 
the  mean  square  error  of  a  reduced  basis  representation  [16] 

M 

xff  =  £  zitVkj,  (9) 

k= 1 

where  M  <  K.  The  most  significant  POD  modes  are  those  corresponding  to  the  largest 
eigenvalues,  since  the  ratio  of  an  eigenvalue  to  the  summation  of  eigenvalues,  A&/  Y,f  At, 
gives  the  percentage  of  the  mean  square  error  unaccounted  for  by  eliminating  the  corre¬ 
sponding  POD  mode  z £  [16]  in  the  reduced  basis  representation.  The  best  stopping  point 
in  the  expansion  (9)  depends,  in  general,  upon  the  application  and  various  algorithms 
have  been  proposed  [19].  For  our  feedback  control/compensator  application,  we  choose 
the  order  M  of  our  reduced  system  to  be  the  maximum  order  such  that  the  associated 
finite  dimensional  control  system  is  both  controllable  and  observable  while  an  acceptable 
level  (99.72%  in  this  case,  see  Section  3  below)  of  data  variability  is  captured. 


2.4  Penalty  Boundary  Formulation  of  the  Control  Problem 


We  use  a  penalty  boundary  formulation  of  the  time-dependent  species  equation  (5)  with¬ 
out  reactions  to  describe  the  transport  of  TMI  in  the  reactor 


dY 

~Bt 


+  v-  vy 

Y (0,  x) 

dY(t ,  x) 

~  dn 
dY(t ,  x) 

~  dn 
dY(t,  x ) 

Bn 


iV  •  (pDVY) 
yo{x) 

i(y(i,  x)  —  u(t))  on  n 
iy  (£,  x)  on  T5 
0  on  T2  U  T3  U  T4  U  T6, 


(10) 


where  Y  is  the  mass  fraction  of  TMI,  W-  denotes  the  outward  normal  derivative,  and 
u(t )  is  the  control.  Under  sufficient  regularity,  one  can  argue  that  solutions  of  (10),  in  the 
limit  as  e  — »  0,  approximate  solutions  for  the  problem  with  conditions  Y (t.  x)  =  u{t)  on 
Tl,  and  Y(t,x)  =  0  at  T5  (see  [4,  10]  for  related  discussions).  For  the  results  presented 
here,  a  value  of  e  =  1  x  10-3  is  used. 


Writing  (10)  in  weak  form  with  test  functions  Wj,  we  obtain 


/ 

Jn 


dY 


~Btw * 


U'.;  dVt  =  — 


(v  ■  vy)  Wj  dil  —  f  DXY  ■  XiUj  dQ 


In  p 


-WjDVY  ■'VpdQ 


WjDY  ds - 


'n,r5 


e  Jri 


WjDu  ds 


(ii) 
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Spatially  dependent  values  for  v.  T,  p,  and  D  are  interpolated  from  the  nodal  values 
obtained  from  the  CFD  simulations. 

As  we  have  explained  above,  two  discretization  formulations  are  applied  to  (11).  The 
first  formulation,  a  finite  element  approximation,  produces  the  full  system  (which  plays 
the  role  of  a  reactor  simulator).  In  this  case,  the  TMI  concentration  in  (11)  is  approxi¬ 
mated  using  standard  finite  element  discretization  with  N  =  453  quadratic  interpolation 
functions  tPi  and  N  nodal  coefficients  yf 


N 

YN(t,  x)  =  .'/f  (/)«'/ (a:0  •  (12) 

i= 1 

Choosing  the  test  functions  to  be  Wj  =  ip j  j  =  1,2,...,  TV,  we  obtain  in  a  standard  way 
the  matrix  equation 

yN(t)=ANyN(t)+BNu(t),  (13) 

where  AN  is  an  N  x  N  matrix,  BN  is  a  vector  of  length  A,  and  the  control  -a  is  a  scalar 
control  function. 

The  second  discretization  formulation  produces  the  reduced  basis  model.  In  this  case, 
we  first  use  the  POD  modes  {zjp}  to  obtain  the  POD  elements 

N 

^k{x)  =  J2zkiipi(x),  k  =  1,  2, ... ,  K,  (14) 

i= 1 

where  the  functions  ipi(x)  are  the  finite  element  quadratic  interpolation  functions.  The 
reactant  mass  fraction  is  approximated  as  a  linear  combination  of  the  POD  basis  elements 
corresponding  to  the  most  significant  POD  modes 

M 

YM(t,x)  =  J2y™(t)Mx)  (15) 

i= 1 

where,  in  this  case,  M  «  K  <<  N.  Application  of  this  approximation  to  (11)  (in  this 
case  we  use  POD  test  functions  Wi  =  i  =  1,2, ,  M)  yields 

yM(t)  =  AMyM(t)  +  BMu(t),  (16) 

where  AM  is  an  M  x  M  matrix,  and  BM  is  a  vector  of  length  M. 
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2.5  Control 


dY 

To  control  the  reduced  order  system  (16),  we  observe  the  flux  of  reactants  —pD—~  to 

on 

the  center  of  the  substrate  ( xp ),  which  is  approximated  as 


M 


qM=-J2pD 


d*k 


k= 1 


dn 


yk(t)  = 


(17) 


We  seek  the  optimal  control  u*  for  (16)  such  that  the  output  qM  tracks  a  signal  qx  (the 
desired  flux  at  xp),  minimizing  the  generalized  performance  index  ( e.g .,  see  [2,  29]) 


V(y0 ,  «(■))  =  /  u'Ru  +  ( q  )'Qiq  +  (q  ~  <tr)'Qi{q  ~  Qt) 


dt  . 


(18) 


Choosing  =  r2i,  R  =  I,  Qi  =  n,  qM  =  (HM)'yM,  (//")'  =  I-  LM(HM)',  and 
Lm  =  Hm ((Hm)' HM)~l ,  we  can  rewrite  this  performance  index  as 


v(yo,  «(■))  =  /  +  (v  ~  yr)'Q(y  ~  v¥)  dt> 


(19) 


where  r i  and  r2  are  design  parameters,  I  is  the  identity  matrix,  Q  =  HM Q2{HMy  + 
HMQi(HMy,  and  y^1  =  LMqT  is  the  desired  state  trajectory.  This  is  a  standard  formu¬ 
lation  for  a  ‘tracking’  control  problem  [2,  29]  and  for  which  a  complete  theory  is  known 
(for  a  summary  and  references  see  Chapter  7  of  [7]).  The  optimal  control  is  given  by 
u*  =  —KMyM  —  gM ,  where  the  gain  is  given  by  (see  p84-85  of  [2],  [29],  or  [7]) 


km  =  R-\BMy  n, 

II  satisfies  the  algebraic  Riccati  equation  (ARE) 

0  =  IL4m  +  (AM)'U  -  UBMR-\BMyU  +  Q, 
and  the  tracking  term  gM  is  given  by 

gM  =R-\BMyb, 

where 

b  =  \(Am  -  fr'/c'1')']  'C 


(20) 

(21) 

(22) 

(23) 


2.6  State  Estimation  in  the  Reduced  Order  Model 


Application  of  the  tracking  control  to  the  reduced  order  model  (16)  yields 


yM  =  AMyM  -  BM KMyM  -  BMg 


M 


qM(t )  =  {. HM)'yM(t ) 


(24) 
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with  KM  and  gM  defined  in  (20)  and  (22),  respectively.  Note  that  this  formulation 
requires  complete  reduced  order  state  feedback  for  the  gain  KM  of  (20).  Since  the  full 
state  of  the  system  (or  even  the  reduced  order  system)  in  the  reactor  can  not  be  observed, 
we  implement  a  state  estimator  or  compensator  design  based  on  observation  (17)  of  the 
flux  at  the  center  of  the  substrate.  This  yields  (see  [29]  or  Chapter  8  of  [7]  and  the 
references  therein)  the  coupled  system 

yM  =  AMyM  -BMKMyf  -  BM gM 
yeM  =  Acyf  +  FM(HM)'yM  -BMgM  , 

where  y ^  is  the  Mxl  vector  approximation  to  the  approximate  state  yM ,  the  compen¬ 
sator  system  operator  Ac  is  given  by 

Ac  =  AM  -  FM(HM)'  -  BmKm ,  (26) 

and  the  compensator  gain  (Kalman  filter)  is  given  by 

Fm  =  EHmV-\  (27) 

The  matrix  £  satisfies  the  dual  ARE  given  by 

AM£  +  £(AM)'  -  T,HMV~l{HM)'T,  +  U  =  0,  (28) 

where  U  is  a  symmetric  positive  semi-definite  matrix  design  parameter  and  V  is  a  sym¬ 
metric  positive  definite  matrix  design  parameter.  We  choose  U  —  I  and  V  =  r3,  where  r3 
is  a  third  parameter  at  our  disposal  in  designing  the  overall  feedback  control/compensator 
system. 


This  implementation  of  the  state  estimator  yields  the  following  closed-loop  system  for 
‘optimal’  control  of  the  reduced  order  model 


/  yM 

V  Ve1 


AM 


M  isM 


-BMK 


FM(HMy  a m _ BmKm _ fm(Hm)'  )  v  yf 


with  the  optimal  control 


u*  =  -KMyf  -gM. 


M 


uM  M  \ 

BM9gM  )  .  (29) 


(30) 


2.7  Control  of  the  Full  System  Using  a  Reduced  Order  State 
Estimator 

Use  of  the  design  given  in  (29)  and  (30)  can  be  expected  to  produce  a  stabilized  and 
generally  efficient  system  for  control  of  the  reduced  order  model  (16).  However,  this  is 
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not  the  issue  of  practical  importance  in  our  efforts.  The  goal,  of  course,  is  to  design  a 
feedback  control/compensator  system  for  (10)  (or  in  weak  form  (11)),  which  we  hope  will 
be  a  good  approximation  (for  proper  choice  of  e)  for  the  actual  physical  dynamics  (l)-(5), 
i.e.,  transient  solutions  of  (5)  with  (l)-(4)  in  steady  state.  Thus,  the  real  measure  of  the 
value  of  the  control/compensator  system  of  Sections  2. 5-2. 6  is  how  well  it  performs  when 
used  in  the  actual  physical  system.  Short  of  applying  the  control/compensator  system 
to  the  physical  reactor  in  experiments,  our  best  assessment  of  its  utility  is  when  applied 
to  the  ‘full  system  (13).  That  is,  we  should  computationally  test  the  reduced  order 
control/compensator  design  based  on  (20)-(23),  (26)-(30)  in  the  system  (13). 

Application  of  the  reduced  order  tracking  control  to  the  full  system  with  the  reduced 
order  state  estimator  yields 

AN  -BnKm  \(  yN  W  -BNgM  \ 

FM(HNy  AM  -  BmKm  -  FM(HMy  H  yM  j  +  l  —BMgM  j  ’  ^  ^ 

where  HN  is  the  observation  vector  for  the  full-dimensional  system  qN  =  (HN)'yN (t), 
and  HM ,  KM ,  gM ,  and  FM  are  as  defined  in  (17),  (20),  (22),  and  (27),  respectively.  The 
optimal  control  is  given  by  (30).  We  report  on  our  simulations  for  system  (31)  with 
different  values  of  the  design  parameter  r\  in  the  results  of  Section  3.  Values  of  r2  and 
t 3  were  varied  and  then  fixed  at  nominal  performance  values  in  the  calculations  reported 
below. 


2.8  Control  Implementation 

The  governing  equations  (10)  are  nondimensionalized  prior  to  the  calculation  of  the  coef¬ 
ficient  matrices  using  a  length  scale  of  1  x  10-2  m,  a  diffusion  coefficient  scale  of  1.15  x  10-3 
m2/s,  and  a  density  scale  of  4.04  x  10-2  kg/m3.  The  POD  modes  and  the  coefficient  ma¬ 
trices  AN ,  BN ,  HN ,  AM ,  Bm ,  and  HM  in  (31)  are  calculated  using  in-house  C  programs. 
All  other  matrix  calculations  are  implemented  in  the  Matlab  computing  environment 
(The  Math  Works  Inc.,  Natick,  MA).  The  optimal  gain  matrix  KM  and  Riccati  solu¬ 
tion  II  are  determined  using  Matlab’s  IqrQ  function  and  the  dynamical  equations  are 
integrated  using  Matlab’s  ode23sQ  solver. 
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3  Results 


We  report  in  sequential  form  results  from  the  series  of  computations  and  simulations 
described  in  detail  in  Section  2. 


3.1  CFD  Simulation  Results 

We  first  report  on  results  of  the  CFD  package  simulations  described  in  Section  2.2. 
Contour  plots  of  the  steady-state  solutions  of  the  temperature  and  the  x-component 
of  the  velocity  (Fig.  2a-b)  give  an  indication  of  the  transport  conditions  in  the  reactor. 
With  the  top  wall  and  inlet  maintained  at  room  temperature,  there  is  a  steep  temperature 
gradient  (Fig.  2a)  upstream  from  the  substrate.  In  the  region  above  the  substrate,  the 
isotherms  run  parallel  to  the  hot  substrate  and  the  opposite  cold  wall.  The  velocity  of 
the  gas  increases  more  than  4-fold  as  it  passes  in  the  vicinity  of  the  hot  substrate  (Fig. 
2b),  while  the  y-component  of  the  velocity  (not  shown)  remains  small  throughout  the 
reactor. 


(a)  Temperature 


(b)  x-component  of  velocity 


(c)  TMI  concentration 


Figure  2:  Contour  plot  of  steady-state  values  for  the  (a)  temperature  in  100°K  steps  with  Tmoa:=750oK 
and  Tmj„=350°K  ,  (b)  x-component  of  the  velocity  in  0.1  m/s  steps  with  vmaa:=0.45  and  v,„m=0.05 
m/s,  and  the  (c)  equilibrium  distribution  of  TMI  2  s  after  the  reactants  entered  the  reactor,  with  contour 
lines  representing  steps  of  0.2  in  the  normalized  TMI  mass  fraction  (Ymoa:=0.9  and  Ymi?l=0.1). 
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A  contour  plot  of  the  equilibrium  distribution  of  TMI  (at  t  =  2  s)  shows  a  steep  gradient 
in  the  TMI  concentration  (Fig.  2c)  at  the  leading  edge  of  the  substrate.  There  is  a 
corresponding  peak  in  the  flux  of  TMI  to  the  upstream  end  of  the  substrate  (not  shown). 
A  plot  of  the  average  flux  of  TMI  to  the  substrate  as  a  function  of  time  (Fig.  3)  shows  that 
an  equilibrium  value  was  reached  approximately  1  s  after  the  reactant  was  introduced 
into  the  reactor. 


Figure  3:  Average  flux  of  TMI  to  the  substrate  as  a  function  of  time. 
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3.2  Construction  of  POD  Modes 


Since  equilibrium  is  reached  after  one  second  (Fig.  3),  only  the  first  sixty  intermediate 
solution  vectors  ( K  =  60),  representing  the  transport  of  TMI  into  the  reactor  from  0- 
1.2  s,  are  used  to  construct  the  POD  modes.  Each  solution  vector  represents  the  TMI 
concentration  at  the  453  nodal  points  and  corresponds  to  a  time  increment  of  0.02-s  in 
the  time  range  from  0  to  1.2  seconds.  A  plot  of  the  captured  variability  CCyLi  \ \) 
as  a  function  of  the  number  of  modes  (M  <  60)  used  (Fig.  4)  shows  that  the  original 
data  is  well-represented  by  7  modes,  or  fewer.  This  strongly  suggests  using  M  <  7  in  our 
reduced  order  model. 


Figure  4:  Total  percent  variability  captured  as  a  function  of  the  number  of  modes. 


The  rank  of  the  controllability  matrix  has  been  found  to  be  a  useful  criterion  (see  the 
discussion  in  [8,  9])  for  determining  the  number  of  modes  to  use  in  control  design  appli¬ 
cations  for  the  reduced  basis  representation  (15).  The  controllability  of  the  linear  system 
(16)  is  determined  from  the  rank  of  the  controllability  matrix  C,  where 

C(Am,  Bm)  =  [Bm  I  AmBm  I  (Am)2Bm  I  ...  I  (Am)m^1Bm],  (32) 

Using  standard  results  from  optimal  control  theory  [14],  the  M-dimensional  linear  system 
(16)  is  controllable  if,  and  only  if,  the  rank  of  the  controllability  matrix  is  equal  to  M. 
Similarly,  the  rank  of  the  observability  matrix  O  indicates  the  observability  of  the  linear 
system  (24),  where 

C(Am,Hm)  =  [Hm  |  (Am)'Hm  |  ((Am)2)'Hm  |  ...  |  ((AM)M_1)'HM].  (33) 
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For  our  system  under  investigation  here,  the  ranks  of  the  controllability  and  observability 
matrices  are  recorded  in  Table  1  as  a  function  of  the  number  of  modes  in  the  reduced  order 
model.  The  Table  demonstrates  that  the  rank  of  the  controllability  does  not  increase  with 
additional  modes  after  addition  of  the  fifth  mode.  Similarly,  the  rank  of  the  observability 
matrix  does  not  increase  with  additional  modes  after  addition  of  the  sixth  mode.  Based  on 
these  results,  we  use  in  our  control/compensator  applications  the  first  five  most  significant 
POD  modes  (M  =  5),  capturing  99.72%  of  the  data  variability,  to  construct  our  reduced 
order  model  (16). 

Table  1:  Rank  of  the  Controllability  and  Observability  Matrices 

Reduce  Order  Model  Rank 
Dimension  M  CO 

1  I  T 

2  2  2 

3  3  3 

4  4  4 

5  5  5 

6  5  6 

7  5  6 
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3.3  Control  of  Full  System  Using  the  Reduced  Order  State  Es¬ 
timator 

We  carried  out  numerous  simulations  in  investigating  performance  of  the  control/compen¬ 
sator  design.  A  set  of  representative  results  are  depicted  here,  where  the  reduced  order 
state  estimator  is  used  to  control  the  full  system  (31)  using  a  non-dimensional  tracking 
value  of  qr  =  0.183  and  design  parameter  values  of  r\  =  100, 180,  or  10000,  r2  =  1  x  10-5, 
and  r3  =  1  x  104.  Initially,  the  solution  vector  for  the  estimated  state  is  given  a  small 
nonzero  value  ( y ^  =  1  x  10-4),  while  the  full  system  solution  is  given  an  initial  value 
of  zero  (yN  =0).  As  a  result  of  the  latter  condition,  the  observed  flux  to  the  substrate 
(Fig.  5)  is  initially  zero  and  this  zero  flux  persists  for  the  time  it  takes  the  reactant 
introduced  at  the  inlet  to  reach  the  substrate.  (Recall  that  the  control  value  represents 
the  reactant  mass  fraction  at  the  inlet.)  As  time  progresses,  the  observed  flux  approaches 
and  oscillates  about  the  target  flux  value.  The  amplitude  of  the  oscillations  increases  as 
the  control  parameter  ri  is  increased,  i.e.,  as  more  weight  is  placed  on  achieving  the 
target  flux.  Eventually  a  steady-state  flux  value  is  attained,  which  exceeds  the  target 
flux  value  for  the  case  of  r\  =  10000  and  undershoots  the  target  flux  value  for  the  case 
of  ri  =  100.  With  a  design  parameter  value  of  r i  =  180,  the  observed  flux  achieves 
the  target  flux  value.  The  steady  state  flux  attained  with  the  control  parameter  value 
=  10000  is  essentially  the  asymptotic  value  for  large  r\,  i.e.,  there  is  essentially  no 
difference  between  the  steady  state  results  for  ri  =  1000  and  r\  =  10000. 


5  10  15  20  25 

TIME  [DIMENSIONLESS] 


Figure  5:  Observed  flux  compared  to  the  desired  flux  (solid  line)  as  a  function  of  non-dimensional  time 
with  different  levels  of  control:  r\  =  100  (dotted  line),  r\  =  180  (dash-dot  line),  r\  =  10000  (dashed 
line) . 
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A  plot  the  norm  of  the  difference  between  the  full  system  solution  and  the  state  esti¬ 
mated  solution  as  a  function  of  time  (Fig.  6)  also  shows  oscillations,  as  initially,  the  state 
estimated  solution  did  not  closely  track  the  full  system  solution.  (Recall  that  the  com¬ 
pensator  or  state  estimator  is  an  asymptotic  estimator  -  see  [2,  7]).  The  greater  the  level 
of  control  (larger  r\)  the  larger  are  the  magnitude  of  the  oscillations.  The  oscillations, 
representing  overshooting  and  undershooting  of  the  adjustments  to  the  state  estimated 
solution,  eventually  (asymptotically)  are  attenuated  and  the  difference  between  the  two 
solutions  reaches  a  steady-state  value. 


Figure  6:  The  norm  of  the  difference  between  the  full  state  and  the  POD  estimated  state  with  different 
levels  of  control:  ri  =  100  (dotted  line),  r\  =  180  (dash-dot  line),  r\  =  10000  (dashed  line). 


A  plot  of  the  control  value  u  as  a  function  of  time  (Fig.  7)  shows  large  initial  control 
values  oscillating  and  decreasing  in  time,  eventually  reaching  steady-state  values.  Both 
the  initial  values  and  the  magnitude  of  the  oscillations  increase  as  ri  increases,  as  more 
weight  is  placed  on  achieving  the  desired  flux  and  less  weight  is  placed  on  the  cost  of 
control.  Physically,  the  control  values  represent  the  mass  fraction  of  TMI  at  the  inlet. 
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TIME  [DIMENSIONLESS] 


Figure  7:  The  control  as  a  function  of  time  with  different  levels  of  control:  r  1  =  100  (dotted  line), 
r\  =  180  (dash-dot  line),  r\  =  10000  (dashed  line). 

4  Conclusion 

We  have  demonstrated  a  proof-of-concept  computational  implementation  of  reduced  or¬ 
der  feedback  control  of  HPCVD  III-V  film  growth  which  can,  in  principle,  be  extended 
to  more  experimentally  relevant  cases  involving  multiple  species  with  reactions.  We  im¬ 
plement  the  control  using  a  reduced  order  state  estimator  based  on  observations  of  the 
flux  of  TMI  at  the  substrate  center,  which  is  precisely  the  type  of  measurements  avail¬ 
able  with  current  PRS  sensing  technology.  The  POD-based  design  method  allows  us  to 
reduce  the  order  of  the  system  by  a  factor  of  90  with  respect  to  a  standard  finite  element 
representation  (from  N  =  453  to  M  =  5),  thus  making  real-time  feedback  control  with 
partial  state  observations  a  feasible  goal  in  HPCVD  reactors  operating  in  steady  state 
flow  regimes  with  pulsed  vapor  reactant  inputs.  This  is,  in  fact,  a  conservative  estimate  of 
the  dimensional  reduction  capability,  since  the  number  of  nodes  used  here  was  minimized 
to  reduce  the  simulation  time  of  the  full  system. 
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